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An efficient scheme is introduced for a fast and smooth convergence to the thermodynamic limit 
with finite size cluster calculations. This is obtained by modifying the energy levels of the non 
interacting Hamiltonian in a way consistent with the corresponding one particle density of states 
in the thermodynamic limit. After this modihcation exact free electron energies are obtained with 
finite size calculations and for particular fillings that satisfy the so called ’’closed shell condition”. 
In this case the ’’sign problem” is particularly mild in the auxiliary field quantum Monte Carlo 
technique and therefore, with this technique, it is possible to obtain converged energies for the 
Hubbard model even for [7 > 0. We provide a strong numerical evidence that phase separation 
occurs in the low doping region and moderate ^ 4t regime of this model. 

PACS numbers: 71.10.Fd, 71.15.-m, 71.30.+h 


After several years of scientific effort, based on ad¬ 
vanced analytical and numerical methods, only very few 
properties of the 2D Hubbard model have been settled. 
The 2D Hubbard model is defined in a square lattice con¬ 
taining a finite number L (Nh) of sites (holes): 

H = K + V = -t 4,aCj,'r + Uj2nlnf (1) 

<ij>,<7 i 

with standard notations, the kinetic energy operator K 
being equal to H for U = 0. In the thermodynamic 
limit, namely for L —>■ oo at given doping S = Nh/L 
fundamental issues such as the existence of a ferromag¬ 
netic phase at large U/t ratio and/or the stability of an 
homogeneous ground state with possible d-wave super¬ 
conducting properties are still highly debated, as several 
approximate numerical techniques lead to controversial 
and often conflicting results. This situation is particu¬ 
larly important right now, since recent progress in the 
realization of fermionic optical lattices could lead to the 
experimental realization of the fermionic Hubbard model. 

Method: We consider the square lattice and generic 
finite clusters that satisfy all the symmetries of the in¬ 
finite system. As well known finite square lattices can 
be defined by two integers n, m such that ri^ -I- = L, 

obtained by supercell translation vectors = (n, m) and 
Ty = (—TO, n). In order to fulfill rotation symmetries, two 
sequences can be defined: 

m = 0 L = n X n The usual sequence (2) 

m = n L = 2v? The 45° degrees tilted sequence 

On the other hand translation symmetries are recovered 
by employing periodic or antiperiodic boundary condi¬ 
tions, the same in both directions and Ty in order to 
preserve rotation and reflection symmetries. 

In the following we would like to consider the most 
useful sequence of clusters for converging as fast as pos¬ 
sible to the thermodynamic limit. A simple technique, 


well known in strongly correlated lattice models^ and in 
realistic calculations^ is to consider the twisted averaged 
boundary conditions (TABC) method. This technique 
allows an exact evaluation of converged thermodynamic 
quantities (energy, density matrix, etc.) in the non in¬ 
teracting U = 0 limit Indeed it has been proven very 
successful to reduce substantially the finite size effects in 
several correlated systems. 

In the following we follow a different approach analo¬ 
gous to TABC in the requirement to remove finite size 
effects in the non interacting limit. However the proposed 
approach can be used more efficiently in combination 
with the auxiliary field quantum Monte Carlo (AFQMC) 
methodi^i^ The latter technique is one of the most power¬ 
ful ones used so far for the study of the Hubbard model, 
as it can project out from a mean-field (Slater determi¬ 
nant) state \MF), the exact ground state of the Hamil¬ 
tonian by the application of the imaginary time propaga¬ 
tion exp(—riJ), for large r. Here we consider variational 
expectation values (Var) on the exp(—r/2i?)|MF) state 
and non variational mixed estimators (no Var) between 
the previous state and exp(—T/2iJ)|^o), IV’o) being the 
ground state Slater determinant for 17 = 0. Both quan¬ 
tities clearly converge to exact values for large t—. 

The sign problem occurs at finite doping and U/t > 0 
but it is particularly mild when i) the U = 0 Hamilto¬ 
nian has a non degenerate ground state, a situation that 
occurs for particular fillings- the closed shell fillings- in 
any finite clusters, ii) the auxiliary field transformation is 
real, as with the proposed approach it is not necessary to 
sample a complex phase. Until now several reliable and 
” numerically exact” calculations have been performed by 
means of this technique on moderately large clusters (i.e. 
up to ~ 100 sites and U/t ^ but it was difficult 

to establish thermodynamically converged results, espe¬ 
cially in the weakly correlated regime. One should also 
mention that, recently, a remarkable progress has been 
made, allowing the complete removal of the time dis- 
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cretization error in the propagato r^i i ^ ^ , a development 
that will not be used here, as we have preferred to use a 
small enough time step Ar, such that this systematic er¬ 
ror is negligible, at least as far as the ground state energy 
is concerned. 

Let us now introduce the method. Consider a finite 
cluster. The kinetic energy can be written in Fourier 
space, by collecting k points related by point symmetries 
to the same energy shell e^: 

^ ~ (3) 

i,<J ki\ek=ei 

where ej, = —2t{coskx + cosky) in 2D. We can assume 
that the energy levels are defined in ascending order. 
Each shell occurs with some multiplicity gi and: 

^ = (4) 

i=l 

where p is the number of different energy levels. In the 
square lattice case we can choose for simplicity clusters 
that do not contain accidental degeneracies {gi < 8), 
namely they are given by the second sequence given in 
Eq.(I2|) with PBC (APBC) for n odd (even). The fun¬ 
damental quantity that we are going to use in order to 
achieve more easily the thermodynamic limit is the den¬ 
sity of states N{E): 


- A (5) 

defined in a way that / N{E)dE = 1. This function 
can be evaluated analytically and/or computed with arbi¬ 
trary accuracy in the thermodynamic limit for any given 
lattice model. We partition the energy bandwidth of the 
lattice (e.g. —At < E < 4t in the square lattice) in inter¬ 
vals ii such that: 


9ilL 


N{E)dE 


( 6 ) 


the following way: 


/ EN{E)dE e, 

= 'A: -= Ljg, f EN{E)dE (7) 

/ N{E)dE -cL 

ei-l 

where the latter equality comes just from the definition 
in Eq.®. In this way the revised kinetic energy K ^ K 
is obtained by replacing Ci with ii in Eq. ([3]) . It is immedi¬ 
ate to show that, when we satisfy the closed shell condi¬ 
tion, i.e. N = J2i<ip 9i^ within these modified boundary 
conditions (MBC) we obtain straightforwardly that the 
ground state energy per site is: 


< K/L >= 2 V —= 2 / EN{E)dE 

^ I 


( 8 ) 


namely the exact energy per site for [/ = 0 at the ther¬ 
modynamic density: 


N/L = 2 j N{E)dE 
eo 

, where the factor two in the above equation takes into 
account the spin components. 

At finite U it is quite simple to show that the above se¬ 
quence of lattices with MBC converges to the exact ther¬ 
modynamic limit because for large L the modification 
of the levels, as compared to the original ones, becomes 
irrelevant. 

In this way we have several advantages and simplihca- 
tions: 

• The Hamiltonian is always real with MBC (a posi¬ 
tive property for the sign problem) 

• The non interacting exact limit is obtained for the 
closed shell fillings. Thus we expect less size effects 
just for those particular densities less affected by 
the sign problem within AFQMC. 


ei-l 

The above equation define all the levels ii by simple in¬ 
duction, because once we know , e.g. we can solve the 
above equation for i = n 1 and we can obtain univo- 
cally in+i- Therefore by setting eg equal to the lowest 
one electron energy {—At in the 2D square lattice), all 
levels ii can be computed and their level spacing is in 
exact correspondence with the density of states. Notice 
that for the particular symmetry of the DOS in bipartite 
lattices N{E) = N{—E), due to particle-hole symmetry, 
it follows that one of the levels is exactly vanishing. 

After the above decomposition, in order to fulfill the 
requirement to have an exact energy for [7 = 0 we can 
modify the energy levels of the kinetic energy Ci ii in 


• The MBC satisfy all the symmetries of the infinite 
systems. For instance when we apply TABC, each 
boundary with a non zero twist generally breaks all 
point spatial group symmetries, maintaining only 
translation symmetry. This may not affect the av¬ 
erage result, but it becomes certainly more diffi¬ 
cult to converge to the exact result, i.e. one needs 
larger projection times in AFQMC due to smaller 
finite size gaps that occur after a small symmetry 
breaking perturbation of the Hamiltonian given by 
a tiny twist of the boundary conditions. 

• Last but not least, MBC are rather trivial to im¬ 
plement in the AFQMC as it is enough to change 
the propagator exp(—ArAT) —>• exp(—ArX). This 
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matrix is never sparse and has to be computed in 
advance within AFQMC for its efficient implemen¬ 
tation. Thus the use of MBC does not lead to any 
overhead in the performances of the algorithm. 
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FIG. 1: U = 0 energy gap for finite clusters with non degen¬ 
erate ground state at the closest filling N/L = 0.9 with 45° 
tilted PBC or APBC. This gap is multiplied by the number 
sites, as it should converge for L —>• oo to g/N{fi), where p is 
the thermodynamic chemical potential at this filling and g is 
the multiplicity of the energy level (g = 8 in 2D and g = 2 
in ID), (a): standard 2D clusters, (b) modified 2D clusters 
according to this work (see text), (c): standard Id clusters, 
(d) Energy per hole for standard clusters with PBC (filled 
symbols) and with MBC (empty symbols). 


As it is shown in Fig.(IT^-c), the most important prob¬ 
lem in the usual sequence of finite clusters is that the 
U = 0 finite size gap behaves erratically when L in¬ 
creases, in spatial dimensionalities D > 1 and away from 
commensurate fillings. This precludes to obtain accurate 
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FIG. 2: Upper panel: Energy per hole in the 2D Hubbard 
model. Finite size effects are rather well behaved by using 
MBC and at small doping the energy per hole approaches 
almost exactly an horizontal line. This implies phase sepa¬ 
ration as discussed in the text. The dashed line is a fit (see 
table) of the data for <5 > 6.5%. Lower panel: energy vs r 
convergence starting with several initial left and right wave 
functions. ”Var” (”no Var”) stands for the (non-)variational 
calculation (see text) with different values of the antiferro¬ 
magnetic parameter Aaf in the mean-field determinant. 


extrapolations to the thermodynamic limit at least for 
moderate U/t and finite dopings (see e.g. FigHJi). As 
it is evident in Fig.([TJi-b) the MBC behave much bet¬ 
ter in this respect. Only few clusters scatter from the 
converged ~ g/N{E) energy level separations, but they 
correspond to atypical lattices when gi ^ 8 at the high¬ 
est occupied or lowest unoccupied free-electron energy 
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levels. Remarkably a reasonably converged value of this 
quantity occurs only after few hundreds sites, which is 
feasible for a numerical approach (see e.g. Fig.[TJl). 

Results: we have carried out a systematic finite size 
scaling study of the energy per hole in the moderate U jt 
regime. As pointed out in the milestone paper by Emery 
and KivelsoiJ^ a minimum at doping (5c in the energy per 
hole, corresponding to a given variational ansatz, implies 
its instability against phase separation because it is pos¬ 
sible to gain energy for (5 < 5c by segregating the holes 
in a hole rich region with the same type of ansatz. In 
an exact calculation, whenever it is possible to carry out 
the thermodynamic limit, clearly we have to obtain a 
constant energy per hole in all the region 5 < 5c, just 
because the compressibility - e.g. the slope in the en¬ 
ergy per hole at 5 —> 0-, cannot be negative in an exact 
calculation. As shown in Fig. [2](a), at C//t = 2 we can 
safely reach the thermodynamic limit with the clusters 
considered, thanks to the very small finite size effects 
introduced by the method proposed in the previous sec¬ 
tion, that, in this case, look considerably smaller than 
standard TABG^. The convergence in imaginary time is 
quite clear also for the most difficult case at small dop¬ 
ing (see Figl^j) and also independent of the initial trial 
function used. In this case we have used a mean held 
state with a non zero antiferromagnetic order parameter 
along the x— spin directioni^ In the more difficult 
cases for larger Ujt^e have also used a d-wave order pa¬ 
rameter in order to converge faster or at least for 

obtaining the lowest possible variational energies com¬ 
patible with a reasonable average sign < s >> 0.05. It 
is clear that this by no means implies the existence of a 
non zero superconducting order parameter in the ground 
state, an issue that will not be discussed in the present 
work. 

At U jt = 2 the achieved flat behavior of the energy 


per hole for 5 < 6.5% clearly indicates the accuracy of 
the AFQMC at this small coupling, that is indeed able to 
determine phase separation just by imaginary time pro¬ 
jection of an homogeneous trial state. At larger coupling 
(see ReflT3l. though we have not been able to reach the 
same cluster size and the same length of the projection 
times, the accuracy in the energy per hole appears ac¬ 
ceptable and allow us to determine the phase separated 
region for U < At (see table). Notice that, in this table, 
the energy gain to have phase separation can be measured 
by the difference of the minimum hole energy ob¬ 

tained for the largest clusters at doping 5 < 5c and the 
hole energy uq extrapolated at 5 = 0 using only dop¬ 
ing values clearly outside the phase separated region. In 
other words uq represents the energy per hole of the uni¬ 
form phase extrapolated at zero doping. This difference 
tto — appears to be very large ~ O.lt at U/t = 2,3 
and less evident for U/t = 4. This is probably the rea¬ 
son why atU/t = 4, there have been several controversial 
claims^i^iSid^. Given this behavior, it is also possible that, 
at larger U/t, the phase separation may be less evident 
and 5c may significantly decrease^^, despite some works 
indicate exactly the opposite effect^di. It is clear that 
at large U/t this important issue remains still open. On 
the other hand, at U/t = 1, we have not obtained evi¬ 
dence of phase separation, because probably we cannot 
reach enough small doping values with the affordable fi¬ 
nite clusters L < 1058. Indeed at this coupling value 
also the antiferromagnetic order parameter ttiaf cannot 
be detected numerically, being extremely small even at 
5 = 0. Since the existence of antiferromagnetism is at 
the basis of the phase separation argumenti^ it is possi¬ 
ble that also 5c can be exponentially small at small U/t, 
as is the case for — exp(— oc \/sJt/U) within the 

Hartree-Fock theory^. 


U/t 

-^min 

5c 

uq 

Ol 

02 

03 

a4 

05 

^max 

1 

-0.484(13) 

0.00(1) 

-0.49185 

0.88825 

1.99156 

-2.46827 

2.20290 

-0.73928 

0.00053 

2 

-0.843(1) 

0.067(5) 

-0.90164 

0.73469 

3.45846 

-4.87765 

4.03989 

-1.27852 

0.00062 

3 

-1.125(1) 

0.105(10) 

1.20909 

0.41515 

4.84506 

-6.43761 

4.74676 

-1.35904 

0.00034 

4 

-1.342(2) 

0.110(15) 

-1.35005 

-0.78626 

9.18032 

-13.03818 

9.60894 

-2.75519 

0.0012 


TABLE I: Estimated energy per hole in the thermodynamic 
limit. The functional form of the fit is a 5*^ order polynomial 
5 

E^{5) = ^ aiS’’ determined by the largest cluster data in 

i=0 

the region 5c < 5 < 1. The rightmost column represents 
the maximum error of the fit for 5 x Eh{S) (the energy per 
site referenced to the undoped case). represents the 

estimated minimum energy per hole for 5^0. Number(s) 
between brackets indicate error bars in the last digit(s). 


Conclusions: We have introduced a technique for con- particularly suited for the AFQMC method. In this way a 
trolling finite size effects in an efficient way, an approach strong numerical evidence is given that phase separation 
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is robust at small dopings and UJt values. With this ap¬ 
proach it is possible to study other possible phasesi^^— , 
with a better control of finite size effects at incommensu¬ 
rate dopings and weak couplings, as well as it is possible 
to export the method to other techniques, that may have 
problem of convergence to the thermodynamic limit es¬ 
pecially at weak couplings. Indeed we have preliminary 
verified that, by choosing boundary conditions that break 
the symmetry of the lattice (e.g. cylindrical), much bet¬ 
ter results (i.e. much closer to the thermodynamic limit) 
can be obtained by correcting the energy levels of the 
[/ = 0 Hamiltonian according to the proposed method. 
Finally we want to remark that the method can be easily 
extended to realistic calculations that do not explicitly 
require a local Hamiltonian^ii^^ as in AFQMC. This can 
be achieved by considering as an input for the correlated 


calculation the band-resolved DOS obtained with an un¬ 
correlated Hamiltonian, such as the Khon-Sham one in 
Density Functional Theory. The same technique as above 
can be used to reduce finite size effects by preserving 
charge neutrality even in presence of the Coulomb long 
range interaction, a property that is difficult to fulfill 
with TABC, if we require that the non interacting limit 
should remain exact with a finite supercell calculation. 
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This supplementary material contains energies of the Hubbard model with the AFQMC method described in the 
paper for various cluster sizes with periodic boundary conditions rotated by 45 degrees. These square lattices have an 
even number of sites L = 2n^ where n is an odd (even) integer with PBC (APBC), and fulfill all rotation symmetries 
of the infinite lattice. In all the tables error bars are between brackets with usual conventions. In all the forthcoming 
results the guiding function is defined by means of the ground state tjjMF of the following mean field Hamiltonian: 


H 


MF 


= K- + 


A 


AF y^(-i) 

R 


k 


( 1 ) 


where K is the modified kinetic energy with MBC and 0 is used when 7 ^ 0 as the non interacting 

chemical potential value, namely at the middle of the non interacting highest energy occupied level tip+i and lowest 
energy unoccupied level Qp, ^0 = 'Rt a^k,a- being the total number of particle operator that is 

fe,(7 

not modified within the MBC approach. R = (x, y) is a lattice point belonging to this lattice, namely R = (x, y) is 
equivalent to (x ±n,y + n). 

The results for the hole energy for U/t = I, 2, 3,4 are summarized in Fig.(IT]). The Trotter discretization time is 
chosen to have a positive systematic error (the values reported are variational upper bound of the energy) in the 
energy per site within 2 x I0“^t, namely At x t = for U/t = 1,2,3,4, respectively. This error, being 

systematic and similar for nearby dopings, is expected to cancel out for the computation of the hole energy at small 
dopings, that, in this case, is affected only by the statistical error. Notice that, thanks to the symmetrized expression 
of the short time propagator: 


exp(—Ari?) ~ exp(—exp(—ArH) exp(—(^) 

the Trotter error in the energy becomes almost negligible if measured after long imaginary time projection, as it is 
scaling as (Ar)"^,— and is also variational in the ”Var” case mentioned in the paper because it corresponds to an energy 
expectation value of a given state. 

For the smallest doping and largest clusters and U/t > 2 values, in order to have converged energies, it is necessary 

to optimize the trial function. Only at U/t = 4 we have found convenient to use a mean field Hamiltonian with a 

2 _ 2 

small d—wave superconducting order parameter. In general A^^g and especially Aaf are useful in the small doping 
region. 


^ At leading order the approximate ground state wave function IV’at) as a function of Ar can be written as IV'At) = + 

At^\iP') + 0(At®) where [ip') is orthogonal to the exact ground state Thus it is simple to show that the approximate 

ground state energy E{At) = does not contain the leading O(Ar^) order in the expansion, and that, with 

WAx IVax ( 

simple inspection, the first non zero contribution is vanishing as Ar'^ 
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L 

Nh 

U/t = 1 

U/t = 2 

U/t = 3 

U/t = 4 

18 

0 

-1.383984(25) 

-1.174357(42) 

-.995837(80) 

-.851727(22) 

18 

8 

-1.316280(18) 

-1.258224(50) 

-1.210676(57) 

-1.171736(17) 

18 

16 

-.404113(4) 

-.402104(12) 

-.400564(15) 

-.399349(4) 

32 

0 

-1.383759(19) 

-1.173957(49) 

-.996540(66) 

-.851421(35) 

32 

8 

-1.427314(18) 

-1.315238(40) 

-1.220670(61) 

-1.141815(28) 

32 

24 

-.802289(8) 

-.791310(14) 

-.782523(24) 

-.775319(17) 

50 

0 

-1.383931(9) 

-1.175175(26) 

-.999904(34) 

-.856989(29) 

50 

8 

-1.432718(9) 

-1.289485(22) 

-1.167422(25) 

-1.065190(50) 

50 

24 

-1.279282(8) 

-1.228975(19) 

-1.187942(20) 

-1.154477(15) 

50 

40 

-.671976(3) 

-.665129(9) 

-.659723(10) 

-.655354(6) 

50 

48 

-.154689(1) 

-.154454(2) 

-.154276(2) 

-.154139(2) 

72 

0 

- 

- 

-1.000386(41) 

-.857926(40) 

72 

8 

- 

- 

-1.123689(54) 

-1.007007(30) 

72 

24 

- 

- 

-1.239880(23) 

-1.180751(37) 

72 

40 

- 

- 

-1.116184(17) 

-1.092080(30) 

72 

48 

- 

- 

-.953056(13) 

-.939974(24) 

72 

64 

- 

- 

-.400496(6) 

-.399212(10) 

98 

0 

-1.383884(17) 

-1.175669(28) 

-1.001029(53) 

-.859353(19) 

98 

8 

-1.417291(113) 

-1.242952(157) 

-1.093107(165) 
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FIG. 1: Energy per hole for various Ujt and cluster sizes. The curved lines are fit to the data with 5*^ order polynomials (see 
main text). The horizontal lines are the lowest hole energies estimated at the smallest available & on the largest clusters. 
























































